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Abstract 

This paper is an introduction to the modelling of viscoelastic fluids, with an em- 
phasis on micro-macro (or multiscale) models. Some elements of mathematical and 
numerical analysis are provided. These notes closely follow the lectures delivered by 
the second author at the Chinese Academy of Science during the Workshop "Stress 
Tensor Effects on Fluid Mechanics", in January 2010. 
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1 Introduction 

Many fluids exhibit a behaviour very different from that of air or water, in particular 
because of some memory effects. Such fluids are called complex or non-Newtonian 
fluids. Examples include blood, muds, fresh concrete, paints, etc... In many cases, 
this peculiar behaviour originates from the presence of some microstructures within 
the fluid the evolutions of which are strongly coupled to the solvent dynamics. In the 
following, we focus on dilute solutions of flexible polymer chains. This is a prototyp- 
ical example of such fluids, both from an experimental and a modelling viewpoint. 
A dilute solution of polymer chains consists of a solvent and polymer chains floating 
therein. These chains are in such a small quantity that direct interactions between 
the chains can be neglected. In the following, a polymer chain is meant to be a long 
linear molecule built as the repetition of an elementary pattern, called the monomer 
(think typically of an alcane molecule CH^ — {CH2)" — CH^). It is observed that 
the rheology of the fluid (that is the way it flows) is very much affected by the poly- 
mer chains, even at a very small concentration. This paper is an introduction to 
microscopic-macroscopic (in short micro-macro) models to describe such fluids. The 
modelling principle is to couple conservation laws on macroscopic quantities (such as 
the velocity or the stress) with some models for the evolutions of microstructures. 
Direct macroscopic approaches will be also discussed. 

We assume that the reader is familiar with standard analytical techniques and 
numerical methods for partial differential equations. Our previous work [57] can be 
seen as a companion paper, more intended to students less knowledgeable on numerical 
and mathematical analysis for partial differential equations. We would also like to 
mention the more general overview fl5' on multiscale methods (not only for polymeric 
fluids, but also for solid materials for example). 

The paper is organized as follows. In Section[2] we present and discuss the building 
blocks of a micro-macro model for a complex fluid, with an emphasis on polymeric 
fluids. This Section contains in particular an introduction to the kinetic modelling of 
polymer chains, and to coarse-graining techniques based on the notion of free energy. 
Section [3] is an overview of the mathematical analysis results that have been obtained 
on models for viscoelastic fluids in general, and micro-macro models in particular. In 
that section, an introduction to entropy techniques for longtime behaviour of partial 
differential equations is provided. Finally, Section [4] is devoted to various problems 
related to the numerical discretization of models for viscoelastic fluids: free energy 
dissipative schemes for macroscopic models, an introduction to a nonlinear approxi- 
mation approach to solve high-dimensional partial differential equations, and various 
techniques for reducing the variance of the numerical results. 

2 Modelling 

The focus of this work is complex fluids whose non-Newtonian behaviour is due to 
some microstructures within the fluid. One example is blood: red blood cells may 
aggregate and form complex structures that affect the rheology of the blood. 
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More precisely, we are interested in cases where: 

1. the niicrostructures are very numerous (per unit volume, say), 

2. the niicrostructures are small and light, 

3. the solvent is Newtonian. 

The first characteristic will be important in the following to express macroscopic quan- 
tities as means over many configurations of the microstructures, using statistical me- 
chanics techniques. The second characteristic will be useful to describe the evolution of 
the microstructures using the Langevin equations, thermalization being modeled using 
a Brownian motion in the evolution equations. Finally, the third characteristic is re- 
quired in order to write the equations at the macroscopic level, assuming that the only 
non-Newtonian contribution to the fluid originates from the microstructures. These 
assumptions do not cover all non-Newtonian fluids with micro-structures, a typical 
counter-example being granular materials, for which other approaches are required. 
However, such a setting is useful for a wide range of materials, a prototypical example 
being dilute solution of polymer chains, which is the main example considered in this 
paper. This setting also provides a convenient test-bed for arguments in mathematical 
analysis and approaches for numerical simulations that may prove useful elsewhere. 

2.1 Experimental observations 

Non-Newtonian fluids are ubiquitous: food industry (mayonnaise, egg white, jellies) ; 
materials industry (plastic or polymeric fluids) ; biology or medicine (blood, synovial 
liquid) ; civil engineering (fresh concrete, paints) ; environment (snow, muds, lava) ; 
cosmetics (shaving cream, toothpaste, nail polish) ; etc... Non-Newtonian fluids may 
exhibit very weird behaviours, such as, for polymeric fluids, the open syphon effect, 
or the rod climbing effect (also called the Weissenberg effect). 

The experimental approach to study Non-Newtonian fluids consists in measuring 
the response of the fluid (the stress) under speciflc constraints (typically homogeneous 
flows), in order to characterize the stress-strain rate relation. A typical experimental 
setting is a simple shear flow obtained in appropriate rheometers, see Figure [T] The 
stress-strain rate relation is obtained by measuring the torque as a function of the 
angular velocity, at a stationary state. The shear stress r is measured for various 
values of the shear rate 7. One can thus discriminate between four typical behaviours: 
Newtonian fluids, for which r = 777, 77 being the shear viscosity ; Shear-thinning 
fluids, for which the effective viscosity T/7 diminishes as 7 increases ; Shear-thickening 
fluids, for which the effective viscosity T/7 increases as 7 decreases ; Yield-stress fluid, 
for which 7 remains until a threshold on the stress t. As explained below, the 
consideration of such simple flows allows to formally classify the non-Newtonian fluids, 
and to propose values for some parameters appearing in macroscopic laws. Of course, 
the challenge is then to assess the validity of these macroscopic laws for more complex 
flows. 

2.2 Multiscale modelling 

Modelling of non-Newtonian fluids starts with the mass and momentum conservation 
equations for incompressible fluids: 



supplied with appropriate boundary and initial conditions. Throughout this article, 
u is the velocity, p the pressure, p is the density of the fluid and cr is the stress. For 
Newtonian fluids, one postulates a linear relation between the stress and the strain 
rate: 




Vp + div(cr), 



(1) 



cr — 



?/(Vm -I- Vu^) 
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Figure 1: Schematic representation of a rheometer. The shear stress is given by r = ^^^2 

where C is the torque and the shear rate is 7 = ^ = ^»ntg — where is the angular 
velocity. 



rj being the viscosity of the fluid. For non-Newtonian fluids, an additional term (the 
so-called extra-stress) is introduced in the constitutive relation 

(T 77(Vm + Vit'^) -1- r. (2) 

The question is now how to define an evolution law on t, as a function of the other 
macroscopic quantities. 

There are basically two approaches for the derivation of such a relation, see Fig- 
ure |2] The first one is to postulate an equation relating r and Vm (using typically 
frame invariance as a guideline) and to fit the parameters to simple experiments, as 
explained above. This leads to two classes of models: (i) integrals models, which pro- 
vide the stress at a given time and a given point as an integral over the past of the 
trajectory which leads to this point: 

T = I mit - t')t(t') dl{t') 

where t(t') is some tensor defined as a function of u (ii) differential models, which 
are partial differential equations involving u and r, a prototypical example being the 
Oldroyd-B model (that we will discuss in more details later on): 

A + M • Vt - Vmt - tVm^^ + t = T]p{Wu + Wu'^), (3) 

where A and r]p are two parameters. This first approach thus yields so-called macro- 
scopic models. They only involve macroscopic quantities (the velocity, the pressure 
and the stress). 

Another route (that we will follow below) is to model the behaviour of the mi- 
crostructures responsible for the non-Newtonian features of the fluid, and to obtain 
the stress r as a function of the configurations of these microstructures. In such an 
approach, additional variables are required (to model the state of the microstructures), 
and one thus obtains higher dimensional systems to simulate. The complete model 
thus couples the Equations ([ij (for mass and momentum conservations on macroscopic 
quantities) to some equations describing the evolution of the microscopic structures. 
This forms a so-called micro-macro (or multiscale) approach. This approach may be 
seen as a way to justify (or find new) macroscopic models, giving a precise substance 
to the parameters that appear in those models. It may also be seen as a modelling 
tool in its own rights, when the micro-macro model does not lead to an equivalent 
macroscopic model, and is discretized as such. 

Multiscale approaches are now very popular for the modelling of materials (flu- 
ids or solids). Such techniques raise several questions, from the modelling viewpoint, 
numerically and theoretically. Theoretically, the main question is the consistency of 
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Figure 2: Classical approaches for the modelling of complex fluids. 



the approach, namely the fact that the models used at various levels are compati- 
ble. More precisely, one would typically like the microscopic models to be refinements 
upon purely macroscopic models. This raises the question of deriving macroscopic 
models from microscopic models (typically using closure approximations, homogeniza- 
tion, etc.). Numerically, many questions are related to the efficiency of the coupling 
procedures (a posteriori error estimators and adaptive refinement techniques to use 
the finest model only when needed, adequate methods to couple models different in 
nature, such as deterministic and stochastic models, etc.). 

From a modelling viewpoint, the delicate question is the discrepancy of the time 
and length scales. Typically, the macroscopic and the microscopic descriptions do 
not share the same time and length scales. The finest time and length scales (at 
the microscopic level) should not be resolved, otherwise the macroscales cannot be 
reached. Two techniques are used. For timescales, one may use local thermodynamic 
equilibrium assumptions to consider that the microscopic configurations reach a sta- 
tionary state (given the macroscopic quantities) within a timescale much smaller than 
the macroscopic timescale. For lengthscales, one may use mean field approximations, 
which consist in replacing the many-particle system at the microscopic level by a one- 
particle problem within an averaged environment. Under these approximations, the 
macroscopic quantities are typically obtained as averages over microscopic quantities 
(locally at the macroscopic time and length scale). Notice that for both time and 
space, these approximations are based on scales separation arguments. Besides, from 
a computational viewpoint, coarse- graining techniques are useful in such approaches to 
reduce the number of degrees of freedom at the finest scale. Typically, coarse-graining 
implies an increase of the microscopic characteristic time and length scales. All these 
techniques will be illustrated below on the modelling of polymeric fluids. 

In the simplest situations, the outcome of the above modelling is a one-way cou- 
pled model (information passing): the microscopic model is simply used to feed the 
macroscopic model, by giving values to some parameters (for example, the viscosity 
of a fluid in the Navier-Stokes equations may be obtained using molecular dynamics). 
But there also exist more complicated situations, where the mean fleld approximation 
to replace the many-body system by a one-body problem at the microscopic level is 
not possible (this is the case of granular materials for example), or where the local 
thermodynamic equilibrium assumption is not valid. In the latter situation, it is not 
possible to assume that the microstructures completely relax to a local equilibrium 
(the micro time scale being of the same order as the macro timescale), and then one 
obtains genuinely coupled systems (information coupling) , as is the case for polymeric 
fluids. 
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2.3 Microscopic models for polymer chains 

Let us momentarily assume that the velocity field u of the solvent is fixed and given. 
The question we address in this section is the modelling of the polymer chain dynamics. 
As explained above, we need a description both sufficiently fine to get the essential 
features of a polymer chain that are of interest to get the correct macroscopic behaviour 
(the correct rheology of the fluid) and sufficiently coarse so that the microscopic model 
is not computationally prohibitively expensive. In the following, we only consider a 
linear chain (without branching). Classical references for the modelling procedure 
described below are [Mj [551 [T^ [TB] . 

2.3.1 Modelling of a chain in the solvent at rest 

We first assume that u — Q. A standard technique to model the evolution of a 
polymer chain is based on molecular dynamics. It amounts to assuming that the 
interactions between the nuclei (these interactions would ideally be obtained using 
electronic structure computations) can be modelled by a potential function, and then 
to writing a dynamics (say, Newton equations of motion) on the nuclei. Let us assume 
that the polymer chain contains N atoms, with positions q — {qi, . . . ,g^) G R"^^, 
interacting through a potential function V . q €z R'^^ i-)- V{q) G JR. In order to avoid 
expensive electronic structure computations, empirical models of the form 

Viqi,...,qN) = ^Vp^iriq^,qJ) + ^triplet (^i, 9^-, ^fc) + •■ • 

i<j i<j<k 

are used for the potential V . For a polymer chain, one would typically have to consider 
pair interactions related to covalent bonds, and interactions involving three or four 
consecutive Carbon atoms along the backbone of the polymer since the potential 
energy depends on angles and dihedral angles along the chain. A standard molecular 
dynamics is then given by the Langevin equations: 

r dQ^^M-^Ptdt, 

\ dPt = -V1/(QJ dt ~ (M-^Pt dt + ^/2C^dWt, ^ ' 

where Pt is the momentum, AI is the mass tensor, C is a friction coefficient and 
/3^^ = kT. Here, Wt denotes a 3A^-dimensional Brownian motion. When C = 0, 
these equations simplify into the Newton equations (Hamiltonian system). The two 
terms depending on C, model the effect of the bath, namely the thermalization due 
to the many-collisions of the polymer chain with the surrounding solvent molecules. 
The term —QM~^Pt dt models a (viscous) friction and the term ^/2CJi~^dW t models 
thermal agitation. These two terms balance each other (the friction rate and the 
diffusion constant are related through the fluctuation dissipation relation) so that the 
canonical Boltzmann-Gibbs probability measure: 

u{dp, dq) = Z^' exp {-p { ^^^^^ + Viq)^ ^ dpdq 

is left invariant by the dynamics. Notice that v{dp,dq) is the product of a Gaussian 
measure on the momentum p with the probability measure 

ii{dq)^Z-^eM-py{(D)dq 

on the positions. Below, we will need the so-called zero mass Langevin dynamics, 
which is the approximation 

= -VF(QJ dt - CdQt + ^/2CjF^dWt 

of Q obtained by letting the mass go to zero. The zero-mass Langevin equation is 
thus a stochastic differential equation on the positions only: 

dQt = -C"' vi/(g,) dt + ^2C~^p-HWt. (5) 
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This dynamics leaves the canonical measure fi invariant. 

This fine description of a polymer chain is not suitable for a multiscale model, 
since it is too high-dimensional. For computational purposes, it is thus desirable to 
derive a coarse-grained model. The assumption underlying the simplest such model, 
the so-called dumbbell model, is that essentially the orientation and the length of the 
polymer chain suffice to describe the influence of the microstructures on the rheology 
of the complex fluid. Let us thus introduce the so-called end-to-end vector 



The function ^ is called a collective variable. To obtain an effective potential, a classical 
procedure is to introduce the so-called free energy Tl{x) which is such that 

$ * (Z-i exp(-/?l/(q)) dq) = exp(-/3n(a;)) dx, (6) 

where ^ * (Z^^ exp(— dq) is the image of the canonical measure by ^. An 
explicit definition is 



Tl{x) = -/3-iln (^j exp(-/?F(q))<5^(g)_,(dq) 



where the measure S^(^q-j_^{dq) (with support {q,£,{q) ~ x}) is defined by the condi- 
tioning formula dq = S^i^q-j_^{dq) dx. The free energy 11 is thus built so that the equi- 
librium thermodynamic properties are exactly conserved, namely, for any (bounded) 
function </? : E,^ — > R 

J ipiaq))Kdq) - J ip{x)eM-mx))dx, (7) 

a property which is equivalent to ^ . For more details on free energy, we refer to [64] . 

For polymeric fiuids, the definition of the potential V to model the polymer chain 
is based upon the Kramers model which describes a freely jointed bead-rod chain. 
Using the coarse-graining procedure described above, one obtains a potential 11 (the 
so-called inverse Langevin potential) which is in practice approximated by: 

• either the FENE potential: 

nFENE(a;) = ^ln(^l-^— ^ 



(iH\x\ 



• or the Hookean potential: 

nHOOK(a;) 

z, 

The parameters H and b depend on the number of segments in the Kramers chain. 
The acronym FENE stands for Finitely Extensible Nonlinear Elastic: the FENE model 
takes into account the finite extensibility of the chain, while the Hookean model is only 
a correct approximation of the inverse Langevin potential around \x\ — 0. 

The dumbbell model thus consists in modelling the chain by two beads linked by a 
spring with potential 11, see Figure |3] Notice that more elaborate coarse-grained mod- 
els have also been introduced, using more than one vector to more precisely describe 
the shape of the polymer chain by a chain of beads linked by springs. The derivation 
we describe here and below apply to such models. 

2.3.2 Modelling of a chain in the moving solvent 

We now assume that the velocity u of the surrounding solvent is still given but non- 
zero. The evolution of the dumbbell (the coarse-grained model of the polymer chain) 
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Figure 3: In the dumbbell model, the polymer chain is modelled by two beads linked by a 
spring. 



consists in writing the zero-mass Langevin on the two beads of the dumbbeU (tlrat is, 
the two ends of the chain). In the following, we use the notation Xl and for the 
positions of the two beads (instead of j and Qjv ti which would be more consistent 
with the latter section) since the beads are thought to be the centers of mass of a few 
monomers, at each end of the polymer chain, rather than the very ends of the chain. 
Then, for i = 1,2, one has 

r dXl - uit, xl) dt = -Q-^F{X\ - X2) dt + y^W^dB], 
\ dXl - u{t, xl) dt = ^Q-^F{Xl - X]) dt + y/2p-^C,-^dBl, 

where B\ and B1 are two independent 3-dimensional Brownian motions. This is the 
zero-mass Langevin equations ([5| written on each bead, using the fact that the friction 
force is now —C^{dX\ — u{t, X\) dt) since the solvent is moving with velocity u, and 
the interaction potential is Tl{xi — X2). In these equations, we have introduced the 
notation 

F{x) = Vn(a;) 

for the so-called entropic force (or mean force) which is the derivative of the free energy 
n. Let us introduce the end-to-end vector 

- - x] 

and the position of the center of mass 

^xl + xl 

* 2 

By linear combination, we obtain 

( dXt = {u{t, X?) - u{t, XD) dt - 2C-^F{Xt) dt + 2^/(FHr^dW], 
\dRt = l in{t, xl) + u{t, X^)) dt + Vp-^-'dWl 

where = ^ [B^ ~ B^j and = [B^ + Bf.) are two independent Brownian 
motions. We now use two approximations: 

• an expansion of the velocity field: 

u(t, XI) ~ u{t, Rt) + Vuit, Rt){X\ - Rt), (8) 

• a smallness assumption on the noise on Rt: 

^/F^Wl « i {u{t, xl) + u{t, X?)) dt, (9) 
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to obtain the following equations: 



dXt = Vu(t, Rt)Xt dt - 2C^F(Xt) dt + ^/Ip^H^dW^ 
dRt = u{t,Rt)dt. 



The Eulerian version of the Lagrangian formulation ( 10 ) is 

dXtix)+u{t, x).VXtix) dt 

= Vu{t, x)Xtix) dt - 2C^F{Xt{x)) dt + ^/IjF^X^dWt 



(10) 



(11) 



In (10 1, the process Xt is (implicitly) indexed by Rq which is the initial position of 
the dumbbell in the solvent. In ( [TT| , Xt{x) depends on the position at time t of the 
dumbbell in the solvent. Notice that the random variable Xt{x) is a function of the 
time t, the position x, and the probability variable lj. 

2.3.3 Discussion of the modelling 

A few remarks are in order concerning the modelling procedure described above. Some 
of these remarks are not specific to the modelling of polymeric fluids, but concern more 
generally molecular dynamics. 

First, concerning the coarse-graining procedure in the solvent at rest (it = 0), we 
have derived the effective potential so that the canonical measure exp{— /3Il{x))dx 
is consistent with the original full canonical measure fj,{dq), at the thermodynamic 
equilibrium (in the sense of ([t])). But then, we have used this potential to derive a 
dynamics on the end-to-end vector: 



dXt = ~2C^\/n{Xt) dt + ^/i|3-\-^dWt. 

It is easy to check that this dynamics indeed leaves invariant the canonical measure 
exp{— f3Il{x))dx (it is in some sense consistent in terms of the thermodynamic equi- 
librium) but a natural question is then the relation between the dynamics (C(Qt))t>o 
and {Xt)t>o- Are the time-correlation functions correctly approximated using the 
coarse-grained dynamics Xt ? Are the transition times between metastable states 
correctly captured ? This is the subject of current researches |62j . 

Second, the coarse-graining procedure is even more questionable for a solvent in 
motion. A flrst question is: what are the correct thermalizing terms in such a situa- 
tion ? In other words, what are the Langevin equations for a particle within a moving 
solvent ? The above approach considers a friction force involving the relative velocities 
—(^{M~^Pt — u{t, Qt)). Can this be justified in a rigorous way ? A second question is 
the relevance of the effective potential 11 originally derived in the solvent at rest for a 
setting where the solvent moves. For example, the image by ^ of the stationary state 
for the original Langevin dynamics on the finest model in a velocity field u ^ is 



certainly not the stationary state obtained for the coarse-grained model ( 11 ). In other 
words, the coarse-graining procedure (using the free energy) is sensitive to velocity of 
the solvent. 

Third, the expansion on the velocity field ([8| introduces the gradient of the velocity 
field in the model while it was originally absent. This leads to some mathematical 
difficulties, since it requires some smoothness on the velocity field. 

Fourth, neglecting the noise on Rt (see (|9|) also has some consequences on the 
mathematical analysis. Keeping the noise term would amount to adding some diffu- 
sion in space on the dumbbell position, which would dramatically simplify the math- 
ematical analysis of the system. 

Finally, we would like to mention that we have described above a suitable model 
for dilute solution of polymers (in particular direct interaction of polymer chains are 
neglected). Similar descriptions have been used to model: 

• rod-like polymers and liquid crystals [Ml 155] , 

• polymer melts [34 1 185 ] . 
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• concentrated suspensions |46] . 

• blood [88] , 

Such kinetic models are thus useful in many contexts. 

2.4 The micro-macro model for polymeric fluids 
2.4.1 Kramers expression for the stress tensor 

In this section, we would like to couple equations ([T| at the macroscopic level to 



equations (111 at the microscopic level. This requires an expression for the extra- 
stress tensor r as a function of the polymer chains configurations. 

Such an expression is standard in the context of molecular dynamics. It carries 
several names: Irving-Kirkwood formula, virial stress or Kramers expression. It writes, 
in our context: 

Tit,x) = np(^ - kTId + ^ {Xt{x) ^ F{Xt{x)))y (12) 

where Up denotes the polymer concentration (which is supposed to be fixed and con- 
stant in this model), E is the expectation (integral with respect to the probability 
variable uj), and (E) denotes a tensor product (for two vectors u and v, v is a 
matrix with elements (uiVj)). This formula can for example be derived starting from 
the definition of the stress in terms of the force exerted on a plane within the fluid, 
see [851. 



2.4.2 The coupled system 

Collecting equations ([TJ, pT| and (12 1, the complete coupled system 
^ m.Vm) = — Vp + t^Au + div(T), 



dt 

dW{u) = 0, 



kTld + 'E{Xt® F{Xt)) 



dXt + uSI^Xt dt = ( \IuXt - -F{Xt) ) dt + 



is obtained. Once non-dimensionalized, it writes: 
' du 



lAkT 



-dWt 



Re 



dt 

div(M) = 0, 



M • Vm = -Vp + (1 - e)Au + div(T), 



We 



(E(Xt(»F(XO)-Id), 



dXt + u ■ W^Xt dt^ [Wu- Xt 



2We 



F{Xt) dt 



/We 



dWt, 



with the following non-dimensional numbers: 



Re=^,We=^,e=^ 



(13) 



where L = y ^ is chosen as the characteristic length scale, \ — is the relaxation 
time of the polymer chains, rjp = UpkTX is the viscosity associated to the polymers, 
and [/ is a characteristic velocity. For the Hookean model, 



F(X) = X and n(X) = 



\xl 

2 ' 



(14) 
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and for the FENE model 
F{X) 



n(X) = -2 In (1 - \X\yb) 



(15) 



Of course, the system of equations (13 1 needs to be supphed with initial conditions 



and boundary conditions. In the following, we denote 2? the physical domain occupied 
by the fluid. 



The micro-macro model (13) can be rewritten equivalently as: 

-Vj3 + (1 - e) Am + div(T), 



Re(^+«.V« 

div(M) = 0, 



Id+ / {X (g) FiX))tP dX , 



(16) 



We 

^ + M . V.^ ^ -divx ( (vmX - ^^FiX) ] ^ 



2We 



where ?/;(<, x, X) dX is the law of the random variable Xt{x). The partial differential 
equation on is called the Fokker-Planck equation. It is derived from the stochastic 



differential equation on X^ in ( 13 1 



2.4.3 Relation to macroscopic models 

It is interesting to notice at this stage that the Hookean dumbbell model is equivalent 
to the Oldroyd-B model (|3|. Indeed, if F{X) = X, it is easy to check that t satisfies: 



+ M.Vr = Vitr + T{Vuf + ^(Vm + (Vtt)^) - — r. 
at We We 



(17) 



Remark 1 The derivative Vt — Vttx — T(Vit)"^ is called the upper convected 

derivative. It is known to be a frame invariant derivative. Other frame invariant 
derivatives exist fl^ [7^, like for example the corotational derivative 

^ +u.Vt -W{u)t -TW{uf, (18) 

where W{u) = . The upper convected derivative is the frame invariant deriva- 

tive which naturally appears in macroscopic models derived from microscopic models. 

There is no known macroscopic equivalent to the FENE model. However, the 
closure approximation [Ml E] E ^ i-\xl\^/b ) — i-E{\xlJ^)/b ^^^^ to the FENE-P 
model, which is another classical macroscopic model (where A = ]E(Xt (E) Xt)): 



e 



We \l-tT{A)/b 



Id 



dA 1 A 1 

— + M • \/A = \7uA + AVm^ , , , H Id. 

dt Wel-tr(A)/6 We 

2.5 Micro-macro models versus macroscopic models 



(19) 



The superiority of the micro-macro approach (system (13)) compared to a classical 
macroscopic approach relies in the quality of the modelling. The approach enables 
numerical exploration of the link between microscopic properties (the polymer chains) 
and macroscopic properties (the rheology of the fluid). This is crucial for deriving more 
predictive models in order to test the behaviour of a known fluid in new regimes, or to 
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MACRO 


MICRO-MACRO 


modelling capabilities 


low 


high 


current utilization 


industry 


laboratories 






discretization 
by Monte Carlo 


discretization 
of Fokker-Planck 


computational cost 


low 


high 


moderate 


computational bottleneck 


HWNP 


variance, HWNP 


dimension, HWNP 



Table 1: Pros and cons for the macro-macro and micro-macro approaches. 



explore the properties of a fluid that has not been synthetized yet. The microscopic 
features of a polymer chain (and more generally of the microstructures within a fluid) 
are indeed typically well known, but the way these microstructures affect the rheology 
of the fluid is very diflicult to predict. Micro-macro modelling is a way to answer this 
question using numerical experiments. 

However, micro-macro approaches remain so far academic tools because of their 
very high computational cost due to the necessary introduction of additional micro- 
scopic variables. Let us also mention that, in our context of polymeric fluids, these 
models share the numerical difficulties of the classical macroscopic approaches. In par- 
ticular, instabilities occur in some geometries when the Weissenberg number becomes 



large (this is the so-called High Weissenberg Number Problem, see Section 4.1). 
These considerations are summarized in Table [l] 

Because of the computational costs mentioned above, it seems natural to try and 
design some numerical methods that combine the macroscopic and the micro-macro 
approaches. The macro-macro model is used where the flow is simple, and the detailed 
micro-macro model is used elsewhere. The idea of adaptive modelling based on a 
posteriori modelling error analysis (see J.T. Oden and K.S. Vemaganti |84j, J.T. Oden 
and S. Prudhomme (83] or M. Braack and A. Ern [24j) has been adapted in this context 
in a preliminary work by A. Ern and T. Lelievre |37| . 



3 Mathematical analysis 
3.1 Generalities 



The main difficulties for the mathematical analysis of the micro-macro systems (13 1 



(or (16)) come (as always) from the nonlinear terms in the equations namely the 
transport terms u ■ Vm, u ■ V dt or u ■ V xipi ^^'^ the coupling terms VuXt dt 
or divx(VMX-0). There exist similar difficulties for the macroscopic models, like the 
Oldroyd-B model. Until recently, only local-in-time existence and uniqueness results 
were available for such models. Recent works |78t [77] by N. Masmoudi (following the 
earlier work |7T] by P.L. Lions and N. Masmoudi) recently changed the landscape. 

The most recent results are based on a priori entropy estimates, and flne results 
on the propagation of compactness in the problem. 

Remark 2 The separation between the coupling terms and the transport terms men- 
tioned above is somehow misleading. Actually, all these terms can be seen as transport 
terms. Let us for example discuss the upper convected derivative which appears in the 



Oldroyd-B model (17). One can check that ify{t,yf{) satisfies: 
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then G{t,y{t,yQ)) = ^y{t,yQ) satisfies 

dG 



dt 



and cr(t,y) = G{t,y)(TQG'^ {t,y) satisfies 
da- 



dt 



u.Vcr — Vucr — (T(Vtt) 



Likewise, for the Fokker-Planck equation in (16 1, one can check that: 
j^{ij{t,yit,y„),Git,yo)X) 

= + «• V,i^ + divx (V,MXi^))(t,y(t,yo),G(t,yo)^)- 

This rewriting of these terms along characteristics is well-known in the literature, 
but it does not seem to be useful as such for theoretical nor numerical purposes (see 
however f61}/). 



3.2 A priori estimates, entropy and longtime behaviour 

A fundamental tool to analyze a model is a so-called a priori estimate which gives 
some bound on the physical quantities. These estimates are useful to prove existence 
results and to analyze the longtime behaviour of the system. In the following, unless 
explicitly stated, we assume homogeneous boundary conditions: tt = on 2?. Notice 
that the a priori estimates we derive and the results on the longtime behaviour we 
obtain below assume the existence of a smooth solution. 



3.2.1 A bad energy estimate 

A first natural energy estimate is obtained by the following procedure. First, multi- 
plying the momentum equation by u and integrating, one obtains: 



Re d 



[ |M|2 + (l-e) / |Vm|2 = - / t: Vm. 
Jt) Jv Jv 



(20) 



Using Ito calculus, the time derivative of the potential energy stored in the springs is: 

(21) 

We recall that 11 is the potential associated to the entropic force F {F — Vll) . Using 
the fact that t — (—Id -I- ]E(Xt (g) F{Xt))) and the the potential 11 is radial (which 
implies that r is a symmetric tensor), one can check that t : Vtt ~ ]E(F(Xf) • VitXt), 
and thus one obtains from (20 1 and (21): 



d 
dt 



Re 

T 



e 
We 



2We 



^E(n(XO)) +(1-6)^|V, 
E(An(Xf)). 



2We^ 



E(|F(XOn 
(22) 



V 



Note that ( 22 1 also writes in terms of the density probability ip (this estimate can 
actually be obtained directly from (16l): 



d 
dt 



Re 
T 



e 
We 



2We 



/ / n^)+(l-6)/ \Vu\^+'f f \F\^^p 

AnV', (23) 
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where is an integral over the microscopic variable X {d is the dimension of the 
ambient space: 

The potential 11 is strongly convex (see ( 14 1 and ( 15 1) so that the right-hand sides 
of (22 1 and (23 1 are positive. This shows that these estimates, although informative 
on a finite time interval, will not be useful to study the longtime behaviour. Since 
no external forces apply on the system, a return to equilibrium is expected: in the 
longtime limit, the velocity should converge to m = and the law of Xt should 
converge to exp{—H). One expects the system to be dissipative and this does 
not seem to translate on ( 22 ) . The question is then how to eliminate the term in the 
right-hand side. For this purpose, we will need to introduce an entropy. 



3.2.2 Entropy estimates: basic facts 

To obtain an energy estimate useful to study the longtime behaviour, we need to 



stochastic counterpart (13) 



consider the micro-macro model ( 16 1 written in terms of the density "0 rather than its 



We first recall the principle of entropy methods ^ on the Fokker-Planck equation 
itself, the velocity field being given. Consider the linear problem: 



Ik 



div; 



-kX 



2We 



Vn(X) 



2We 



(24) 



where tp is only a function of the time variable t and of the microscopic variable X, 
and K is a given tensor, which is the gradient of the velocity field. A now classical 
approach to study the longtime behaviour of such an equation consists in introducing 
the relative entropy 

±_ 

."000 



ff(^l^oo) 



In 



(25) 



of -ip with respect to a stationary state ipoo of ( 24 1 . Notice that H > and that H = 
if and only ii ijj = tpao- Let us recall the Csiszar-Kullhack inequality 



which shows that the entropy ( 25 ) provides a control on the L^-norm of ■;/; — -00 
the fact that -^oo is a stationary solution to ( 24 ) , one can check that 



rfg(V'|0oo) 

dt 



1 



2We 



Vln 



0. 



(26) 



Using 



(27) 



The right-hand side is non positive, and the entropy H thus decreases in time. 

In order to obtain the convergence of to with a specific rate of convergence (and 
thus the convergence of ijj to ipoo using (26)), one then needs a functional inequality 
on the stationary measure V^oo, called a logarithmic Sobolev inequality. It is said that 
■000 satisfies a logarithmic Sobolev inequality if there exists a constant p > such that 



In 



1 

< — 
- 2p 



Vln 



(28) 



for all non-negative function with integral 1. Inserting (28) in (27) and using the 



Gronwall Lemma, one obtains the exponential convergence: > 0, 



H{i:{t, OIV-oo) < H{^{0, OlV'oo) exp (- 



We 



We are left with the question of proving a logarithmic Sobolev inequality ( 28 1 for 
^poo- This is a delicate question (especially when no analytical expression for tp, 



IS 



available) and a very active subject of research. A related issue is the identification 
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of the spectrum of a Fokker-Planck operator with a drift term that is not in gradient 
form. 

Let us recall some standard criteria to obtain a logarithmic Sobolev inequality. 
There currently exist three basic tools to prove that a probability measure satisfies a 
logarithmic Sobolev inequality. The first one is the Bakry- Emery criterion: 

Proposition 1 Let 11 he an a-convex function, that is a function such that for all 
vector X, X'^V^UX > a\X\^ (wh ere V 11 is the Hessian of II). Then, the density 
probability ■i/'oo ct exp(— 11) satisfies a logarithmic Sobolev inequality with constant 
p > a. 

The second one is a perturbative result due to HoUey and Stroock: 

Proposition 2 Let 11 6e a function such that the probability density ipao oc exp(— 11) 
satisfies a logarithmic Sobolev inequality with constant p. Let II be a bounded function, 
and let us introduce the probability density oc exp(— II+II). Then tp^ also satisfies 
a logarithmic Sobolev inequality with constant p > pexp(— osc(n)) where osc(n) = 
supll — inf n. 

The third tool will be less useful in our specific context but is crucial in other fields. It 
relates to tensorization of probability measures. If {i'i)i<i<i are probability measures 
that satisfy logarithmic Sobolev inequalities respectively with constant pi, then vi (g) 
. . . <^ vj also satisfies a logarithmic Sobolev inequality, with constant min(pi, . . . , pj). 
This result for a collection of / independent random variables admits generalizations 
when some correlations are introduced, see |87l HTl 155] . 



Going back to the Fokker-Planck equation (24 1, one can check (SHI that: 

• If /« is zero or skew-symmetric, ip converges exponentially fast to ipoo oc exp(— 11) 
(since 11 is a-convex). 

• If K is symmetric, an analytical expression for the stationary state is again avail- 
able: ipoo cc cxp{—Il + WeX^ kX) , and thus one obtains exponential convergence 
for the FENE model for any values of the Weissenberg number We (using the fact 
that in the FENE model, X leaves in a bounded domain), and for the Hookean 

model under the assumption that the eigenvalues of k are strictly smaller than 

1 

2We- 

• For a general tensor k and using HoUey-Stroock criterion, exponential decay 
to equilibrium is obtained if there exists a stationary state V'oo that satisfies 
II ln(-!/)oo exp(n))||icxj < oo. This can be proved for the FENE model under the 
assumption that the symmetric part of k is sufficiently small (see ([33| below) . 



3.2.3 Entropy estimates for the coupled micro-macro system: u 
on dV 



We now return to the coupled problem ( 16 1, and first assume for simplicity that u — 



on the boundary dV (no external forcing) . In this natural and expected station- 

ary state is Uao = and -000 oc exp(— H). Let us reconsider the energy estimate (23 1, 
and introduce the entropy instead of the energy stored in the springs. This consists in 

replacing in ( |23p the term Hip by In (^^^^ "0 (which is, up to an additive 

constant, Hip + iplmp). Then, a straightforward calculation yields: 



d_ 
dt 



2We 



Vvln 



F{X)-VuXiP, 



■D o/R" 



(29) 
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instead of pTj) . By linear combination of ( 20 ) and ( 29 1 , one obtains (compare 



with (231) 



di 



Re 

T 

(1 



e 
We 



In 



iVwl 



e 

2 



V'oc 



2We" Jv JR"* 



Vxln 



V'oo 



(30) 



= 0. 



The energy ^ 



V 

Re f |„,|2 



Wc /d /r-^ IIt/i is not a decreasing function (see (p3|), but the 
We /d /r-^ ('^) ^ indeed decrease in time (see (30)). 



free energy ^ Jjy\u 

The dissipative nature of the system is obvious on (|30 1 



Using the logarithmic Sobolev inequality satisfied by -^oo (which holds since 11 is 
a-convex), and the Poincare inequality 



/ |Mp < Cpi(2?) / \Wu\ 
Jv Jv 



(31) 



one obtains the convergence of {u, tp) to the equilibrium (1X007 "^oo) = (0, Z''^ exp(— 11)) 
at exponential rate. 

Two remarks are in order. First, in order to write the free energy of the system, 



we have considered the Fokker-Planck version ( 16 1 of the micro-macro system and not 



the stochastic formulation ( 13 1. It is unclear how to obtain a similar estimate only in 



terms of the random variable Xt{x) (which would be useful for analyzing the stability 
of numerical schemes based on the stochastic formulation (13 1, for example). Second, 
on the linear problem (24 1, it is actually possible to analyze the longtime behaviour 



using many entropy functionals of the form 



ir(^|V'oo)= h 
Jr" 



where /i is a regular scalar convex function such that h{l) = ^'(1) = see [J. Clas- 
sical examples include h{x) = a;ln(x) — a; + 1 (classical "physical" entropy, related to 
logarithmic Sobolev inequalities) and h{x) — (a;— 1)^ (x^-distance, related to Poincare 
inequalities). On the nonlinear problem (16 I, h{x) 



X ln(x) 



1 which yields the 



entropy ( 25 1 introduced above is seemingly the only choice that leads to the dissipa- 
tive structure (30 1. This is a generic observation: for many nonlinear Fokker-Planck 



equations (in the sense of MacKean), the correct entropy seems to be the "physical" 
relative entropy ( |25[ ). This may be related to the fact that this entropy naturally 
passes to the limit when considering discretizations using interacting particle systems 
(extensivity property). 



3.2.4 Entropy estimates for the coupled micro-macro system: u ^ 
on dV 

A natural question is to generalize the result of the previous section to a system 
with non-zero forcing. This falls into the general question of convergence to non- 
equilibrium steady states. The question is known to be much more difficult than the 
study of equilibrium steady states (namely steady-stated without external forcing, at 
the thermodynamic equilibrium). 

In our context, we consider a forced system through the boundary: it 7^ on 
&D. In such a situation, a stationary state (mqcV'oo) is not the equilibrium state 



(0, exp(— n)), but is simply defined as a stationary state to (16 1 satisfying the 
boundary conditions. For simplicity and without loss of generality, let us assume in 
this section that: Re = 1/2, We = 1 and e = 1/2. 
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Mimicking the reasoning of the previous section we obtain, instead of (30 1, 



dt \2jj, 



'"'if I* 



iVtt 



Vxln 



= - / M-VttooM- / / M • Va;(ln-0oo)V' 
I'D JV JR'' 



(Vx(ln^oo) + Vn(X)) • VmXV, 



(32) 



where u{t, x) = u{t, x) — Urxi{x) and ijj{t, x, X) — ijj{t, x, X) — ^00(3^7 Compared 
to the situation m = on dfl, there are two new difficulties: (i) new terms in the 



right-hand side of (32 1 appear and have to be controlled; (ii) one needs to prove that 
■0OO satisfies a logarithmic Sobolev inequality. 

For the FENE model, it is actually possible to obtain exponential decay to equi- 
librium. A typical situation where a proof can be performed is the following (we refer 
to [50 for more details). We assume that the boundary conditions are compatible 
with a homogeneous stationary state: Uoo{x) = kx (where /« is a trace free matrix). 
Then, V'oo is assumed not to depend on the space variable x. Thus, the second term 



in the right-hand side of (32) is zero. The first term can be controlled by the left-hand 



side if /t is sufficiently small, using a Poincare inequality on u. To have an estimate 
of the last term, as well as to prove a logarithmic Sobolev inequality for ip^o (using 
Proposition |2| , one needs to control |lVx lii(V'oo exp(n))||ioo. This can be achieved 
under an additional assumption on k. More precisely, one can prove that if 

< 1/2, 

where = ^(k + k^), then there exists a unique stationary solution ipoo which 
satisfies: VX e 8(0, Vb), 



V In 



exp(-n(X)) 



- 2k'X 



< 



1 - 2k^| 



(33) 



where [., .] is the commutation bracket: [ Using this result, 

one may prove that u converges exponentially fast to Urx, in the L^-norm, and the 



entropy / / \ ; 

assumption on k: 



4, 



tp converges exponentially fast to 0, under an additional 



M^&^ exp(46Af) + C'pi(2?)|k:"| < 1 



where Cpi(V) is the constant in the Poincare inequality (31 1 and M = 2|k'' | + -j^i 



3.2.5 Entropy estimates for macroscopic models 

In this section, we are interested in macroscopic models, such as the Oldroyd-B model, 
or the FENE-P model. The aim is to analyze the longtime behaviour of such models, 
using ideas developed in the similar study for micro-macro models. In the following, 
we assume that m = on dT>, but some generalizations in the spirit of the latter 
section are possible. 

We have seen in Section p3. 2. 3| that the analysis of the longtime behaviour of micro- 
macro models relies on the free energy ^ jitp + Jx> /r<i (^^) ^ rather than 

the energy ^ \u\^ + ^ Jj^Jj^a^'ip- Since we mentioned that some macroscopic 
models (such as the Oldroyd-B model for example) are equivalent to some micro- 
macro model, it is thus natural to ask how the results obtained for micro-macro 
models translate to macroscopic models. We refer to [651 [TTl 1551 llOll H7] for more 
details. 
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We first consider the Oldroyd-B model ([s]), whiclr is equivalent to the Hookean 
dumbbell model. To simplify the presentation, let us rewrite the model in terms of 
the so-called conformation tensor 



A 



We 



r + Id 



(34) 



which writes in the micro-macro formulation as the positive symmetric tensor: A 
E{Xt (81 Xt). The Oldroyd-B model rewrites in terms of {u,p, A): 



Re 



du 

m + 

dW{u) = 
dA 

1h 



u Vu] = -Vp (1 - e)Au + —div(A), 
' We 



M • VA = VuA + AVm^ (A- Id). 

We 



(35) 



To mimic the derivation of ( |22[) , we. on the one hand, multiply by u and integrate 
over V the equation on u in (35) (which yields ([20|): 



Re d 



/^H' + ,i-„/jv„|--4/^. 



A : Vm 



(36) 



and, on the other hand, take the trace of the equation on A and integrate it over T): 



— / tr A = 2 / A : Vm 



V 



1 

We 



tr A + —d\D\, 
X, We ' 



(37) 



where d denotes the dimension and the Lebesgue measure of T). By linear combi- 
nation of (36 1 and (37), we obtain: 



d /Re 



dt 



V 



2We 



J trA^+(l-e) J \Wu\ 



2We^ 



trA = 



2We- 



;d\V\. (38) 



This equality is equivalent to (22) for Hookean dumbbells (n(X) = |Xp/2). It thus 
shares the same drawbacks as (22) for the study of the longtime behaviour, related to 
the positive right-hand-side in ( 38 ) . 

To write a free energy estimate of the type ( 30 1 , one first establishes the additional 
identity: 

d 



dt 



f trlnA=-^- / tr(A-^-ld). 
Jv We 7p ' 



(39) 



This identity is obtained from the equation on A in ( 35 ) , by contraction with A ^ 
and integration over V, using the Jacobi identity: 



d_ 
Ft 



u V]a] : A-^ 



d_ 
di 



M • V (tr In A). 



(40) 



Combining (|36|, (|37| and (39), we obtain: 

,9 e 



d /Re 
dt 



V 



il-e) j \Vu\ 
Jv 



2We 
2 I 



/ tr (A -In A -Id) 
Jv / 

I tr(A + A"^ -2Id) =0. 
Jv 



(41) 



2We 



Equation (41) is equivalent to (30). In particular, using that, for any symmetric 
positive matrix M, 

< tr(M-lnM-Id) < tr (M + M"^ - 2Id), 



and using the Poincare inequality ( |3l[ ), one obtains the exponential convergence of 
(m,A) to the equilibrium states (moo, Aqo) — (0,ld). 
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This result on the Oldroyd-B model may be employed as a guideline to derive 
suitable energy estimates useful for the study of the longtime behaviour of other 
macroscopic models, even if they are not known to be equivalent to a micro-macro 



model of the form (13 1. For the FENE-P model (19) for example, one obtains 



e f f tr(A) 2d ( a-^\\ 

'^m? JvW- tr " 1 - tr {A)/b + V~ 

(42) 

Using next the Poincare inequality (31 1 and that for any symmetric positive matrix 
M (with size d x d), 

< - ln(det(M)) - 6 In (1 - tr (M) /b) + (b + d) In 



b + d 



, tr(M) 2d 1. 

< ^ htrfM"^ 

- \(l-tr(M)/6)2 l-tr(M)/6 ^ 



one obtains that the free energy 
Re 



X ''''' + 2We X A) - Mn (1 - tr {A)/b)) 

decreases exponentially fast to 0. 

Besides the study of the longtime behaviour, the interest of such estimates is 
twofold. They can be used to propose new stability criteria for numerical schemes 



(see the notion of free energy dissipative schemes introduced in Section 4.2 below). 
These estimates can also be used to establish new existence results (see Section 3.3.2[ ). 



3.3 Existence results 

We would like to now present some results concerning existence and uniqueness of mod- 



els for viscoelastic fluids. Notice beforehand that systems like (35 1 modelling purely 
macroscopically non-Newtonian fluids typically include the Navier-Stokes equations, 
with the additional term divA in the right-hand side. The equation on A is essentially 
a transport equation and, formally, A has at best the regularity of Vit. The term 
divA in the right-hand side of the momentum equation is thus unlikely to bring more 
regularity on u. It is then expected that the study of these coupled systems contains 
at least the well-known difficulties of the Navier-Stokes equations. Recall that for 
the three-dimensional Navier-Stokes equations, it is known that global-in-time weak 
solutions exist but the regularity, and thus the uniqueness, of such solutions for ap- 
propriate data is only known locally in time. We will see below that basically, the 
picture is similar for some models for viscoelastic fluids. 

As mentioned in Section [3J] the difficulties in the mathematical analysis of such 
models are related to the transport terms (in addition to u ■ Vm, we now have u ■ VXt 
and u ■ Vijj), the nonlinear terms either coming from the coupling between variables, 
or inherently contained in the equations defining the stress (due to the non-linear 
entropic force F in micro-macro models, for example). 

3.3.1 Existence results for partial differential equation formulations: 
macroscopic models 



The state-of-the-art of the mathematical knowledge for the Oldroyd-B model (35 1 
is existence of local-in-time strong solutions. Example of results for such macro- 
scopic models are contributions by M. Renardy in [S3] (local-in-time existence for 
an abstract model that covers the specific Oldroyd-B case), by C. Guillope and 
J.C. Saut |43l Hi] (existence results for less regular solutions for non-zero viscosity 
of the solvent), by E. Fernandez-Cara, F. Guillen and R.R. Ortega ^39) and references 
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therein (local-in-time well-posedness in Sobolev spaces), by F.-H. Lin. C. Liu and 
P.W. Zhang [68_ (local-in-tinie existence and uniqueness results and global-in-time 
existence and uniqueness results for small data). 

The first global-in-time existence result is due to P.-L. Lions and N. Masmoudi |71| . 
where an Oldroyd-like model is considered, with the corotational derivative ([T8| rather 
than the upper convected derivative on the stress tensor (see Remark [l] above) . The 
proof relies on the fact that for the corotational derivative, some estimates can be 
obtained on the stress, independently of u. 

In a recent work [77 , N. Masmoudi shows global-in-time existence of weak solutions 



to the FENE-P model (19 1, and other nonlinear macroscopic models (Giesekus model 



and Phan-Thien- Tanner model) . The proof is based on the entropy estimate ( 42 1 , and 
a very fine analysis of how compactness propagates in the system. 

3.3.2 Existence results for partial differential equation formulations: 
micro-macro models 



We now turn to the micro-macro models in their Fokker-Planck formulation ( 16 ). The 
analysis of such micro-macro models for polymeric fluids has begun with an early work 
by M. Renardy [M] , 

The difficulties present for the purely macroscopic models and discussed above are 
also present mutatis mutandis in the multiscale models. 

Several existence and uniqueness results for the coupled system involving the 
Fokker-Planck equation have been obtained. For some local existence and uniqueness 
results, we refer to M. Renardy [94j, T. Li, H. Zhang and P.W. Zhang [66^, H. Zhang 
and P.W. Zhang |103j . N. Masmoudi [75] . The latter author has also shown a global 
in time existence result for initial data close to equilibrium (see also F.-H. Lin, C. Liu 
and P.W. Zhang [69 for a prior, more restricted result). 

Global existence results have also been obtained for some closely related problems: 

• Existence results for a regularized version: In a series of paper [7l|8l|9], J.W. Bar- 
rett and E. Siili obtain global existence results using space regularized versions. 
For example, the velocity u in the Fokker Planck equation is replaced by a 
smoothed velocity, and the same smoothing operator is used on the stress tensor 
Tp in the right-hand side of the momentum equations. Alternatively, a small 
diffusion term (with respect to the space variable) is introduced in the Fokker- 
Planck equation. See also L. Zhang, H. Zhang and P.W. Zhang |104| . 

• Existence results with a corotational derivative: In J.W. Barrett, C. Schwab and 
E. Siili \ix |8i (again with some regularizations) and P.-L. Lions and N. Mas- 
moudi [m |79] (without any regularizations), the authors obtain global-in-time 
existence results replacing Vit in the Fokker-Planck equation by '^"^J^" . More 
precisely, in [72] , a global-in-time existence result of weak solutions is obtained 
in dimension 2 and 3, while in [73], it is proved that in dimension 2, there exists 
a unique global-in-time strong solution. A related recent result by F.-H. Lin, 
P. Zhang and Z. Zhang is [70]. 

The most striking existence result has been obtained by N. Masmoudi recently, 
who was able to prove the global existence of weak solutions for the FENE model, 
see [78]. 

We would like also to mention the related works |281 [2^ [3D] by P. Constantin, 
C. Fefferman, N. Masmoudi and E.S. Titi, on existence results for coupled Navier- 
Stokes Fokker-Planck micro-macro models. 

3.3.3 Existence results for stochastic differential equation formula- 
tions 



We now turn to system (13 1, which is quite difficult to study in full generality. For 



preliminary arguments, a standard simplification is to consider simple shear fiows (see 



20 



in particular in M. Laso and H.C. Ottinger J.C. Bonvin and M. Picasso |18| . 
C. Guillope and J.C. Saut [44J, W. E, T. Li and P.W. Zhang [35 ). In such a case, 
the velocity is in the a;-direction, and only depends on the y-variable. Such flows are 
typically studied in rheometers, see Figure |4] 





velocity profile 

Figure 4: Schematic representation of a rheometer. On an infinitesimal angular portion, 
seen from above, the flow is a simple shear flow (Couette flow) confined between two planes 
with velocity profile {u{t,y),0). 



In this simple setting, the system (13 1 writes: 



^ du , , , d u , , dr 



dPtiv) 



^{t,y)Q,{y)dt-^Jp{Xt{y))dt 
dQtiy) = -^FQ{Xdy)) dt + -^dWf, 



(43) 



/We 



2We 



/We 



where t is the off-diagonal component of the stress r, {Pt{y),Qt{y)) are the two 
components of the stochastic process Xt{y), and {Fp,Fq) the two components of the 
force F. 

For Hookean dumbbells, the global-in-time existence and uniqueness is easily ob- 
tained, since Qt is known independently of {u, Pt), so that the system becomes linear, 
see the work by B. Jourdain, C. Le Bris and T. Lelievre \52\ . 

For the FENE force, two new difficulties have to be addressed: first, the stochastic 
differential equation contains an explosive drift term and second, even in a shear 
flow, the coupling term VuXt is genuinely nonlinear. We refer to B. Jourdain and 
T. Lelievre |5T] and to B. Jourdain, C. Le Bris and T. Lelievre |52| for the first studies 
in this setting and in particular the proof of a local-in-time existence and uniqueness 
result. A surprising result is that if b (the maximum extensibility) is not large enough, 
the entropic force is not sufficient to prevent the dumbbell to go out of the centered 
ball with radius Vb. 

For a similar result in a more general setting (3-dimensional flow) and forces with 
polynomial growth, we refer to W. E, T. Li and P.W. Zhang [33]. The authors prove 
a local-in-time existence and uniqueness result in high Sobolev spaces. We also refer 
to A. Bonito, Ph. Clement and M. Picasso |17j for existence results for Hookean 
dumbbells, neglecting the advection terms. 
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More generally, it is important to notice that, when the velocity field is not suffi- 
ciently regular, and similarly to the situation seen for the Fokker-Planck setting above, 
it is difficult to give a sense to the transport term in the stochastic differential equation 
(which is actually a Stochastic Partial Differential Equation). We refer to C. Le Bris 
and P.-L. Lions [591 150] for works in this direction. 



4 Numerical methods and numerical analysis 

Most of the numerical methods employed in practice to simulate models for viscoelastic 
fluids are based upon a finite element discretization in space (see however R. Owens 
and T. Phillips [H^ for spectral methods) and a finite difference discretization in 
time (usually Euler schemes), with a decoupled computation of {u,p) and t. More 
precisely, at each timestep, the equation for {u,p) is first solved, given the current 
stress tensor t. This allows to update the velocity. Next, the equation for t (possibly 
the underlying stochastic model) is solved, and the stress is updated. 



4.1 Generalities 

Of course, the difficulties raised by the discretization of the models always, 
reminiscent of the difficulties of the mathematical analysis. Here again, as mentioned 
above, the treatment of the multiscale problem necessarily requires a good knowledge 
of the treatment of the purely macroscopic model. An overview of the numerical 
difficulties encountered when simulating purely macroscopic models for non Newtonian 
fluids may be found in R. Keunings [ST, F.P.T. Baaijens [S], and R. Owens and 
T. Phillips [89 . 

Let us mention in particular: 

• An inf-sup condition is needed between the discretization space for r and that 
for u (in the limit e — f). This requires either appropriate discretization spaces 
(satisfying some inf-sup conditions, see |76| ) or stabilization methods (see |42) 

m)- 

• The discretization of the advection terms needs to be adequately performed. 
Typically, one may use again stabilization methods [40j or a numerical charac- 
teristic method |20l If02l [50] . These difficulties are prominent for high Reynolds 
number (which is not practically relevant in the context of viscoelastic fluid 
simulations) or for high Weissenberg number (which is relevant). 

• The discretization of the nonlinear term raises difficulties. In particular, treating 
in the equation on the stress the velocity explicitly and the stress implicitly leads 
to an ill-posed problem if the Weissenberg number becomes large. 

We have mentioned that two of these difficulties are prominent for large Weis- 
senberg number. It indeed appears that numerical methods become unstable in this 
latter regime. This is the so-called High Weissenberg Number Problem (HWNP) al- 



ready mentioned in Section 2.5 Many works are related to the HWNP (we refer for 
example to R. Owens and T. Phillips [13 Chapter 7]). The HWNP is certainly not 
only related to the discretization scheme. It has indeed been observed that for some 
geometries, the critical Weissenberg number (above which the scheme is unstable) 
decreases with the mesh step size (see R. Keunings [S^), which could indicate a loss 
of regularity for the continuous solution itself (see D. Sandri [HS])- It is still an open 
problem to precisely characterize the HWNP, and to distinguish between instability 
coming from the model itself, or its discretization. For the theoretical study of the 
limit We — )■ oo, we refer to M. Renardy [95l Chapter 6]. 

In the following, we mainly focus on the discretization techniques for micro-macro 



models. There are two approaches: discretizing the Fokker-Planck version (f6l, or 



discretizing the stochastic version (13 1. Let us start with a section on macroscopic 



models and free energy dissipative schemes. 
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4.2 Discretization of macroscopic models: free energy dissipa- 
tive schemes 

In this section, we would like to discuss one peculiar aspect of the discretization 
of macroscopic models, such as the Oldroyd-B model or the FENE-P model. The 
question we would like to address is the following: what are the necessary ingredients 



such that a numerical scheme for the Oldroyd-B model ( 35 1 has similar dissipative 
properties in terms of free energy as the continuous version, see (41 1. This would 
ensure longtime stability of the numerical scheme. 

This question has interesting practical consequences since it is often the case that 
the stationary state is obtained as longtime computation on the time-dependent prob- 
lem. Besides, as already mentioned above, the discretization of models for viscoelastic 
fluids is known to be difficult due to some stability problems when the Weissenberg 
number becomes large, at least in some specific geometries. These unstable behaviours 
may have two origins: the model itself, or the discretization method. In the first case, 
the unstable behaviour would be related to some singularities of the continuous so- 
lution which appear in finite time. For the Oldroyd-B model for example, it seems 
that the model itself is sometimes questionable, see |92l IMl E]. This may be seen 
as a counterpart of the fact that the Hookean dumbbell can extend to infinity, a 
highly unphysical fact. Other cases are unclear. A numerical scheme recently pro- 



posed in |38l I3H] and based on a reformulation of ( 35 1 in terms of the logarithm of the 

tensor A seems to stay stable for Weissenberg numbers larger than the usual thresh- 
old observed for standard numerical schemes. In this context, one aim of the study 
of the free energy dissipative structure of the numerical scheme is to help to distin- 
guish between stability problems due to the numerical scheme, and stability problems 
originated form the model itself. 

Let us introduce the standard variational formulation used for a finite element 



discretization of the Oldroyd-B model (35 1: find {u,p,A) such that for all {v,q,B 



[ ( Re I ^ -t- u ■ Vit I ■ V + (1 — e)Vu : Vv — pdivv + =^A : Vv + gdivit 
Jv \ \ dt J We 

(44) 

a A \ 1 \ 

— -t- M • VA (WuA + AVm^) :B+ (A - Id) : B =0. 

dt J We / 



To obtain the free energy estimate (41 1, one in principle just needs to take (v, q, B) = 
(^u,p, 2i|ro(Id — j4~^)) test functions. This is obviously not possible for any arbitrary 
finite element space. 

A numerical schemes for a which an estimate similar to ( [4l] ) holds at the discrete 
level is the following: using Scott- Vogelius finite elements ((Pd)^ x IPd-i.disc in di- 
mension d = 2 or d = 3) for the velocity-pressure, Pq finite elements for the tensor 
A, and a discretization of the advection terms by the characteristic method, or a 
discontinuous Galerkin method. As an example, in the latter case, the variational 
formulation is: find {u'l+\p'l+\ A';^+^) € (Fd)^ x lPd-l,d^sc x (Pq)^ such that for all 
[v, q, B) e (Fd)2 X Pd-iMsc X (Po)': 

+ ul ■ ^ul+^ ]-v- pl+^d:iYV + qdivM^'+i 




(45) 



{Vul+^Al+' + Al+\Vul+^f) : B j 
\ul-n\{Al+'\:B+ 

where {Kk)i<k<NK denote the mesh cells, {Ej)i<j<NE the internal edges of the mesh, 



Ne 

E 
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n is the normal on the edge. the downstream value of (/) (for the vector field 
u]l) and (Al+^j = (^1+^)+ - {Al+^y is the jmnp of Al+^ through the edge. The 
following stability result then holds: 

Theorem 1 For a given mesh, there exists a constant co > such that for all timestep 



sizes 6t < Co, there exists a unique solution to (45 1, with A^^'^^ a positive definite 



matrix. Moreover, this solution satisfies the following free energy estimate: 

(46) 



where 



StJ^ (^(1 - e)|V<+ip + ^tr (A^i + {Al+Y' 27rf)) < 



We refer to [23' for more details on free energy dissipative schemes. The main conclu- 
sions therein are: 

(i) The discretization method of the advection terms plays a crucial role in the 
proofs. A numerical characteristic method or a discontinuous Galerkin method 
are adequate discretization techniques. 

(ii) The discretization space for the stress should contain piecewise constant func- 
tions. This is related to the fact that Pq is a stable finite element space under 
the application A i-> 2We(Id — (which is the test function we would like to 
choose to recover the free energy estimate). 

(iii) It is possible to obtain a free energy dissipative scheme without limitation on the 
timestep size using the logarithmic formulation of the stress, as recently intro- 
duced in [38l [48] (whereas a CFL condition is required for standard formulations 
in terms of A, see Theorem [l]). Moreover, any solution to the logarithmic for- 
mulation satisfies a free energy estimate (and this is not true for the standard 
formulation in A). This is maybe related to the fact that these schemes seem 
to be more stable than standard ones, but definite conclusions in this direction 
require further studies. 

4.3 Discretization of the Fokker-Planck version: solving high- 
dimensional partial differential equations 



Spectral methods are typically used to discretize the Fokker Planck equation in (16 1 
(see A. Lozinski [73] or J.K.C. Suen, Y.L. Joo and R.C. Armstrong [HZI)- It is not 
easy to find a suitable variational formulation of the Fokker Planck equation, and 
an appropriate discretization that satisfy the natural constraints on the probability 
density ip (namely non negativity, and normalization). We refer to C. Chauviere and 
A. Lozinski [27] [71] for appropriate discretizations in the FENE case. One major 
difficulty in the discretization of Fokker-Planck equations is when the configurational 
space is high-dimensional. 

The main difficulty with this approach comes from the high-dimensionality of the 
function 'ijj. For the dumbbell model in dimension 3, tp is a, function of 7 scalar vari- 
ables. When the polymer chain is modelled by a chain of N beads linked by springs, 
the Fokker-Planck equation is a parabolic equation posed on a 3A^-dimensional do- 
main, and ip is a, function of 4 -I- 3A^ scalar variables. Some numerical methods have 
been developed to discretize such high dimensional problems. The idea is to use an 
appropriate Galerkin basis, whose dimension does not explode when dimension grows. 
We refer to P. Delaunay, A. Lozinski and R.G. Owens |32], T. von Petersdorff and 
C. Schwab [lOOj . H.-J. Bungartz and M. Griebel [55] for the sparse-tensor product 
approach, to L. Machiels, Y. Maday, and A.T. Patera [7S] for the reduced basis ap- 
proach and to A. Ammar, B. Mokdad, F. Chinesta and R. Keunings [U |2] for a 
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method coupling a sparse-tensor product discretization with a reduced approximation 
basis approach. 

We more precisely describe the latter method, which is an interesting nonlinear 
approximation method. For simplicity, we concentrate on the Poisson equation in 
dimension 2, but the generalization to any symmetric problem on a high-dimensional 
cylindrical domain is straightforward. 

Consider the problem: 

Find gei/i(f^) such that I "^'^"i ™ ^' (47) 
' [ g = on ail, ^ ' 

where / e L^(i7), 17 = x with fi^ C H and ilj, C R two bounded domains. The 
bottom line of the algorithm proposed in [1] [2] is to look for the solution as a sum 
of tensor products g = I]fe>i '^fc ® (where rfc (8) Sk{x,y) = rk{x)sk{y)) computing 
iteratively the terms in the sum: set /o /, and iterate on n > 1, 

1. Find r„ e Hl{Vlx) and s„ e Hl{Q,y) such that: for all functions (r, s) G Hl{Q.x) x 
Hl{Vty) 

V(r„ (8) s„) • V(r„ (g) s -I- r (8) s„) = / /„_i(r„ s -1- r (8 s„). (48) 
n Jo, 

2. Set /„ = / + A^"^irfc(8sfe. 

3. If ||/n||ff-i(j2) > £, then go to iteration n -|- 1, otherwise stop. 

Notice that this method can indeed be generalized to high dimensional Poisson prob- 



lems since the problem ( 48 I in dimension N consists in solving a (non- linear) system 
of N equations. This has to be compared with a classical Galerkin approach, which 
would lead to a linear problem with dimension that exponentially scales in N . This 
method reportedly gives reliable results in high dimension, see 2-^ and [81^ J2^ for 
applications in the context of uncertainty quantification. 

A fundamental remark to understand the algorithm is that the variational prob- 



lem ( 48 1 is the Euler equation associated to the minimization problem: Find r„ S 
Hliyix) and s„ G Hl{VLy) such that, 

(r„,s„) = arg ^ min \\ I W {r ® s)\'' ^ ( f°_^{r ® s)] . (49) 



'(r,s)effi(n,)xHi(o„) V2 



Then, as shown in [58 , the proposed algorithm falls within a class of so-called Greedy 
Algorithms introduced in nonlinear approximation theory, see [101 1311 1531 198| . Using 
results from ||33|, it is possible to prove the convergence of the algorithm with rate 
0(l/y/n). For convergence results for nonlinear problems (minimization of strongly 
convex functionals) , we refer to [26 . 

The convergence of such algorithms for non-symmetric problems is, to the best of 
our knowledge, an open question. 

4.4 Discretization of the stochastic version 



As mentioned in the previous section, discretizing the Fokker-Planck version ( 16 1 of 
the micro-macro models is certainly not the most natural approach especially when 
the microscopic model involves numerous degrees of freedom (a chain of springs rather 
than a dumbbell to model the polymer chain, for example). 

4.4.1 The CONNFFESSIT method 



The standard approach to discretize ( 13 1 couples a finite element discretization and 
a Monte Carlo technique. In this context, this approach is called CONNFFESSIT 
for Calculation Of Non- Newtonian Flow: Finite Elements and Stochastic Simulation 
Technique (see M. Laso and H.C. Ottinger [56]). 
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To discretize the expectation in (13 1, a Monte Carlo method is employed: at each 



macroscopic point x (i.e. at each node of the mesh once the problem is discretized) 
many replicas (or realizations) [X^' )i<k<K of the stochastic process Xt are sim- 
ulated, driven by independent Brownian motions {W'l)k>i, and the stress tensor is 
obtained as an empirical mean over these processes: 




F(X^")-Idj. (50) 

This approximate stress is then inserted in the mass and momentum conservation 
laws, to update the velocity and the pressure. 

One important feature of such a discretization is that, at the discrete level, all 
the unknowns {u,p,t) become random variables. The consequence is that variance 
is typically the bottleneck for the accuracy of the method. In particular, variance 
reduction methods are relevant and important. 



4.4.2 Convergence of the CONNFFESSIT method 



For simplicity, we concentrate on a simple shear flow (43). In this case, and for 



Hookean dumbbells, the CONNFFESSIT method writes as follows: at time t„ = n5t 
{6t being the fixed timestep size) , for a given velocity field G Vh and dumbbell con- 
figurations {Pl'",Q'ii")i<k<K e dyVh, compute £ Vh and {P^"''^^ ,Qh'^~^^)i<k<K e 
dyVfi solution to: 



r Re 



e 1 
We^ 



k,n+l 



P 



j^k,n+l 



Here [W^ ^Wf^") i<k<K is a collection of independent two-dimensional Brownian 
motions, Vh is a finite element space (think of PI finite elements, namely continuous 
piecewise affine functions) and dyVh denotes the space of derivatives with respect to 
y of functions in Vh (think of PO functions, namely piecewise constant functions). 
For such a problem, a typical result of numerical analysis first proved in B. Jourdain, 
C. Le Bris and T. Lehevre [SIJ (see also W. E, T. Li and P.W. Zhang |3Sj) is the 
following error estimate: 



P 



k=l 

•k,n 



1 -k.n 



2We 



2We 



/We V 



u{t„ 



E(p*„gtj--^E^^'"^?'^'" 

fe=l 



1 



where Ay is the size of the mesh associated to the finite element space V/j. The 
precise result actually requires a cut-off procedure, which is applied with very small 
probability in the limit At — )■ or X — )■ c», see |52] for the details. 

The main difficulties for the proof originate from the following facts: 

• The velocity u\ is a random variable. The energy estimate at the continuous 
level cannot be directly translated into an energy estimate at the discrete level 

due to some correlation between random variables. The stability of the scheme 

k 

is thus not trivial to obtain, and actually requires an almost sure bound on Qh 
which is obtained through a cut-off procedure. 



The end-to-end vectors {Ph 



4h 



)i<k<K are coupled random variables (even 



though the driving Brownian motions {Wf'^,W^''^)i<k<K are independent). 
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For an extension of these results to a more general geometry and a discretization by a 
finite difference scheme, we refer to T. Li and P.W. Zhang [671 . A convergence result 
in space and time may be found in A. Bonito, Ph. Clement and M. Picasso 1 116] . 

4.4.3 Variance reduction for the CONNFFESSIT method 

We would like to now discuss two variance reduction methods for the CONNFFESSIT 
discretization. 

The control variate method The first method is the control variate method. 
Let us start by a general presentation of the method. 

Assume that, in full generality, one wants to compute the expectation ]E(Z) of a 
random variable Z using the Monte Carlo method. The idea of the control variate 
method is to write 

E(Z) = E(Z - oY) + aE(y) 

(where the so-called control variate, is a random variable with a known average, 
and a is a deterministic parameter to be fixed) and to approximate the expectation 
E(Z — aY) using a Monte Carlo method. For the approach to be efficient, we need 
the control variate Y to be such that 

Var(Z-ar) < Var(Z). 

The optimal value of a can be explicitly computed. It is a* = '^°var(r)^'* that 
Var(Z - a*Y) = Var(Z)(l - p{Z, Y)^) with 

KZ,Y)= e[-M]. 

v/Var(r)Var(Z) 

This shows that the optimal case is y)| = 1 (i.e., up to an additive constant 
Y = Z — E(Z)) but this requires to know E{Z) which is precisely what we want to 
compute. The worst case is p{Z, Y) = (that is Y and Z are uncorrelated) in which 
case the control variate does not help in reducing the variance. The practical question 
is then how to build a control variate Y sufficiently correlated to Z, but such that 
E(y) is explicitly known (or sufficiently easy to compute) ? 

In the context of micro-macro models, the Monte Carlo method is used to approx- 



imate the stress r (see (50 1 ) and one can build control variates using two companion 
systems: (i) the same problem at equilibrium or (ii) the same problem with a force 
for which a macroscopic equivalent is known (Hookean force, or the FENE-P closure 
approximation, for example). For the FENE model for example, one may write: 

-f E(Xt®F(Xt)) , 

with either: 

• (control variate at equilibrium) FiX) = F{X) = i_\x\^/b ^"^^ '^■^t ~^ ^ 

yx,dt = ~^^F{X,)dt + ^JWt. 

• (Hookean dumbbell as control variate) F{X) = X and dXt + u ■ WXt dt = 



{VuX,^^^F{X,)) dt+^JW,. 

In both cases, E (^Xt (X) F{Xt)j can be computed using deterministic approaches: 

it is either the stress at equilibrium (thus an analytically known constant), or the 
solution to the Oldroyd-B partial differential equation. Of course, for the method to be 
effective (namely for the control variate to be indeed correlated to the original random 
variable), the Brownian motion driving Xt needs to be the same as the Brownian 
motion driving Xt- 
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Dependency of the Brownian motion upon the space variable Another 
variance reduction technique is very specific to micro-macro models such as (13 1. It 



consists in using the space correlation of the driving Brownian motion as a way to 
control the variance of the result. 

We again consider Hookean dumbbells in a shear flow (see ([43|): 



^ du , , 



e] — {t,y) 



r{t,y) 

dPtiv) 
dQt = 



e 

We 

du 

dy 
1 



2We 



^{Pt{y)Qt) 
{t,y)Qt dt- 
Qt dt 



dr 
dy 



it,y), 



2We 



Ptiy)dt 



(51) 



--dwy, 



One important remark (valid in any geometry and for any force) is that the macro- 
scopic variables (m,p, r) do not depend on the correlation in space of the Brownian 
motions. This owes to the fact that t is defined as a pointwise average. At the discrete 
level, however, the variance of the results does depend on this correlation, which thus 
may be seen as a adjustable numerical parameter that may help reducing the variance 
of the results. 

We study this problem in on the simple system (511. In particular we justify 



theoretically the numerical observations |45l [T5] : using on different mesh elements 
mutually independent Brownian motions increases the variance on the velocity and 
decreases the variance on the stress compared to using a Brownian motion constant in 
space. Loosely speaking, the fact the variance on the velocity is large with mutually 
independent Brownian motions is related to the fact that a derivative with respect to 
space of the Brownian motion appears through the term ^ in the right hand side of 
the equation on u. In some formal sense, this is a derivative of a white noise. The 
fact that the variance on the stress is smaller with mutually independent Brownian 
motions is related to the fact that the stress appears to be an integral in space of the 
Brownian motions, and a sum of independent random variables as a smaller variance 
than a sum of the same random variable (namely Var ^X][=i < Var 
where G* are i.i.d.). The main conclusions of are: 

• The variance of the results comes from an interplay between the space discretized 
operators and the covariance of the Brownian motion in space. 

• The minimum of the variance of u is obtained for a Brownian constant in space. 

• The minimum of the variance of r is not obtained with mutually independent 
Brownian motions. One can further reduce the variance using a Brownian motion 
Wt multiplied alternatively by +1 or —1 from one cell to another (in the spirit 
of variance reduction through antithetic variates). 



4.4.4 A stochastic reduced basis approach for micro-macro prob- 
lems 

We conclude this article summarizing a recent work |22| on a reduced basis approach for 
Monte Carlo computations. The motivation is the following: in the CONNFFESSIT 
method, the expected values: 

TEiXt(g>F{Xt)) 
where Xt satisfies the stochastic differential equation: 

dXt + u ■ \7^Xt dt^ (\/u- Xt —F(Xt)] dt + -^^dWt 

V 2We ^ *7 /We 

have to be computed through a Monte Carlo procedure for many values of the ve- 
locity gradient. The principle is to use this "many query" context to build a variance 
reduction method. 
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The two building blocks are: the control variate method already explained above 
in Section [4.4.3| and the reduced basis method, which we now outline. 

The reduced basis method The reduced basis method [211 |2I] has been pro- 
posed to efficiently solve partial differential equations, typically of the form — div(yl(A)Vu) = 
/, when the solution u{X) has to be repeatedly computed for many values of the pa- 
rameter A e A. 

The method consists: 

• in an offline stage, in building a reduced basis 

Xjv = span(u(Ai), . . . , u{\n)) 

for appropriately chosen values of the parameter A, and 

• in an online stage, in looking for a variational solution to the original problem 
in Xn, using a Galerkin method. 

In the offline stage, the parameter values Ai are computed using a greedy algorithm: 
choose a finite set Atrial C A, and Ai G Atrial- Then, for n > 0, set 

A„+i e arg sup A„(A) 

AeAtHal 

where A„(A) is an a posteriori estimator of the error made on the output, when 
approximating u(A) by Un{\) € span(u(Ai), . . . , u(A„)). 

We purposely omit to mention other aspects, and in particular the fact that, to be 
efficient, the method requires some orthogonalization procedure of the reduced basis. 
This approach has proven to be useful in many contexts. A satisfactory accuracy is 
typically achieved for N of the order of 10. In summary, the main tools used in the 
reduced basis methodology are: 

• A two-stage offline- online strategy (which is efficient either in a many-query 
context, or for real-time computing); 

• the idea to use solutions at given values of the parameter to build a reduced basis 
on which solutions for other values are expanded; 

• A procedure to select the best linear combination on a given reduced basis; 

• An a posteriori estimator used online and offline to evaluate the error. 

• A greedy algorithm to select offline the best values of the parameter among a 
trial sample to build the reduced basis; 

We now describe a method, based upon the same ideas, to build a control variate to 
reduce the variance for Monte Carlo estimations by empirical means of parametrized 
random variables. 

A reduced basis method for Monte Carlo computations Consider the 
problem of approximating, for many values of A, ]E(Z^) by 

M 
i=l 

where {Z^,Y^) are i.i.d. We use the empirical variance 

M 
i=l 

as an a posteriori estimator of the error. As mentioned above, an optimal, unfortu- 
nately inaccessible, control variate is = — ]E(2''^). 
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We thus propose the following two-stage algorithm. In the offline stage, optimal 
control variates are computed, using a large number of samples: 

-'"large .^j^ 

In the online stage, we use 

n=l 

where (a*) — arginf(c,^') Var ^Z"^ — ^^^^ a„F'^"^ , as a control variate. The op- 
timal (a*) are solution to a least square problem. Then, the expected values are 
approximated online using empirical means over a small number of samples (typically 
-Wiargc — lOOMsmaii)- Finally, the parameters A„ in the offline stage are computed 
using a greedy algorithm, in order to minimize the variance. 

In [22], we apply this technique to the FENE model, ]E(Z'^) being the stress tensor 
and the parameters A being the components of the gradient of the velocity field. One 
representative numerical result is displayed on Figure |5] As seen on this Figure, the 
conclusions are: 

• Using the control variate method with a 20-dimensional reduced basis of control 
variates, the variance is approximately divided by a factor of 10* in the mean for 
large test samples of parameters, at least in the applications we experimented. 
As a consequence, the reduced-basis approaches allows to approximately divide 
the online computation time by a factor of 10^, while maintaining the confidence 
intervals for the output expectation at the same value than without reduced 
basis. 

• The greedy algorithm seems to efficiently select the values of the parameter A 
since in the online stage, on a large set A of values for the parameter A, the 
reduction on the variance is basically of the same order as the one observed in 
the offline stage, on the trial sample Atrial- 




H) 2 4 6 S W 12 14 16 18 K) 2 4 6 8 10 12 14 16 18 20 



Figure 5: Maximum, median and minimum of the variance over the samples, in the offline 
stage, over Atnai (left) and in the online stage, over A (right), as a function of the size of 
the reduced basis. 

There exists another version of the algorithm specialized to the case when Z'^ is a 
function of some solution to a parametrized stochastic differential equation. We refer 
to [22] for more details. 
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